Electric displacement as the fundamental variable in electronic-structure calculations 
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Finite-field calculations in periodic insulators are technically and conceptually challenging, due to 
fundamental problems in defining polarization in extended solids. While significant progress has been 
made recently with the establishment of techniques to fix the electric field £ or the macroscopic 
polarization P in first-principles calculations, both methods lack the ease of use and conceptual 
clarity of standard zero-field calculations. Here we develop a new formalism in which the electric 
displacement D, rather than £ or P, is the fundamental electrical variable. Fixing D has the intuitive 
interpretation of imposing open-circuit electrical boundary conditions, which is particularly useful 
in studying ferroelectric systems. Furthermore, the analogy to open-circuit capacitors suggests an 
appealing reformulation in terms of free charges and potentials, which dramatically simplifies the 
treatment of stresses and strains. Using PbTi03 as an example, we show that our technique allows 
full control over the electrical variables within the density functional formalism. 
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The development of the modern theory of polariza- 
tion [l[ has fueled exciting progress in the theory of the 
ferroelectric state. Many properties that could previ- 
ously be inferred only at a very qualitative level can now 
be computed with quantum-mechanical accuracy within 
first-principles density-functional theory. Early ab-initio 
studies focused on bulk ferroelectric crystals, elucidating 
the delicate balance between covalency and electrostat- 
ics that gives rise to ferroelectricity. Over time, these 
methods were extended to treat the effects of external 
parameters such as strains or electric fields 0, Of 
particular note is the recent introduction by Dieguez and 
Vanderbilt Q of a method for performing calculations at 
fixed macroscopic polarization P. The ability to compute 
crystal properties from first principles as a function of P 
provides an an intuitive and appealing link to Landau- 
Devonshire and related semiempirical theories in which 
P serves as an order parameter. 

Despite its obvious appeal, however, the constrained-P 
method has found limited practical application to date. 
One reason for this is that the procedure to enforce a con- 
stant P during the electronic self-consistency cycle is rel- 
atively involved; this hampers the study of complex het- 
erostructures with large supercells, where computational 
efficiency is crucial. There are also physical reasons. In 
particular, fixing P does not correspond to experimen- 
tally realizable electrical boundary conditions (Fig. [I}. 
Moreover, in an inhomogeneous heterostructure, the lo- 
cal polarization can vary from one layer to another, and 
its average is therefore best regarded as a derived, not 
a fundamental, quantity. In the following we show that 
considering D as the fundamental electrical variable over- 
comes these physical limitations, and that constraining 
D rather than P leads to a simpler implementation. 

Formalism. We consider a periodic insulating crystal 
defined by three primitive translation vectors a^, with f2 





FIG. 1: Electrical boundary conditions within differ- 
ent methods, a) The fixed-£ method corresponds to adopt- 
ing closed-circuit boundary conditions with a constant applied 
bias V. b) Constraining D corresponds to a capacitor in open- 
circuit conditions with a fixed value of the free charge Q on 
the plates, c) constraining P does not correspond to a clear 
experimental set-up. 
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the unit cell volume, and we introduce the new functional 

n 

8tt 1 

U(D, v) depends directly on an external vector parameter 
D, and indirectly on the internal (ionic and electronic) 
coordinates v through the Kohn-Sham energy Eks and 
the Berry-phase polarization P [l[ . (For the moment we 
fix the lattice vectors; strains will be discussed shortly.) 
The minimum of U at fixed D is given by the stationary 
point where all the gradients with respect to v vanish, 
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Comparing with the fixed-£ approach of Ref. @, Q 
which the electric enthalpy T is given by 
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F(£,v)=E KS (v)-ne--p(v) , 
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provided that we set £ = D 47rP. Wc thus discover 
that D = £ + 4irP is the macroscopic electric displace- 
ment field. The functional in Eq. (Q]) takes the form 
U = £ks + (O/871") £ 2 , which is the correct expression for 
the internal energy of a periodic crystal when a uniform 
external field is present (details are given in Supplemen- 
tary Section 2.4). Eq. (QJ thus provides a framework 
for finding the minimum of the internal energy U (D) 
with respect to all internal degrees of freedom at speci- 
fied electric displacement D. This is the essence of our 
constrained-D method. 

As a consequence of Eq. (j4|) , the method is analogous to 
a standard finite-£ -field calculation [2, [3(. In particular, 
the Hellmann-Feynman forces are computed in the same 
way. The only difference is that the value of £, instead of 
being kept constant, is updated at every iteration until 
the target value of D is obtained at the end of the self- 
consistency cycle (or ionic relaxation). This implies that 
the implementation and use of the constrained-D method 
in an existing fmite-£ -field code is immediate; in our case 
it required the modification of two lines of code only. 

The effect of constraining D, rather than £, essentially 
corresponds to the imposition of longitudinal, rather than 
transverse, electrical boundary conditions. For exam- 
ple, as we shall see below, the phonon frequencies ob- 
tained from the force-constant matrix computed at fixed 
D are the longitudinal optical (LO) ones, while the usual 
approach yields instead the transverse optical (TO) fre- 
quencies. Furthermore, the longitudinal electrical bound- 
ary conditions are appropriate to the physical realization 
of an open- circuit capacitor with fixed free charge on 
the plates, while the usual approach applies to a closed- 
circuit one with a fixed voltage across the plates (Fig.Q}. 

Stress tensor. This analogy with an open-circuit ca- 
pacitor suggests an intuitive strategy for deriving the 
stress tensor, a quantity that plays a central role in piezo- 
electric materials. In particular, the electrode of an iso- 
lated open-circuit capacitor cannot exchange free charge 
with the environment. This suggests that the flux of the 
vector field D through the three independent facets of 
the primitive unit cell should remain constant under an 
applied strain. These fluxes are x & ■ D = f2 bj ■ D, 
where the are duals (aj • bj = <£y ) differing by a factor 
of 27T from the conventional reciprocal lattice vectors. We 
then rewrite the functionals in terms of the "internal" or 
"reduced" variables di = (f2/47r)bj • D. It is also use- 
ful to define the reduced polarization pi = b^ • P and 
the "dual" reduced electric field £j = a$ • £ . Additional 
details are provided in the Supplementary Section 4. 

By Gauss's law, di = —Qi, where the Qi are the free 
charges per surface unit cell located on the cell face nor- 
mal to hi. With these definitions, the internal energy 
can be rewritten as 

U{{d}) = E KS + ^ " Pi) 9ij ( d j ~ Pj) (5) 



where we have introduced the metric tensor gij = a, ■ a.j 
We then define the fixed- {d} stress tensor as 
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where 77^ is the strain tensor. By a Hellmann-Feynman 
argument (sec Supplementary Section 4.3) the total 
derivative in Eq. ([6]) can be replaced by a partial deriva- 
tive. Using dgij/dr/nv — 2aj tt o 7 - v , we find 



r KS 1 „Max 1 „aua 



where 0™ is the standard zero-field expression, 
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is the Maxwell stress tensor (which originates from the 
derivative acting on and fi _1 ), and 



dpi 



(9) 



is the "augmented" part. If the internal variables v are 
chosen as reduced atomic coordinates and plane- wave co- 
efficients in a norm-conserving pseudopotential context, 
neither the ionic nor the Berry-phase component of pt 
has any explicit dependence on strain, and u™ s vanishes. 
The name thus refers to the fact that <r^ g is nonzero only 
in ultrasoft pseudopotential [B| and projector augmented- 
wave @ contexts. 

We note that, as a consequence of fixing the reduced 
variables di rather than the Cartesian D, the proper 
treatment of piezoelectric effects 0, 0] is automatically 
enforced. This formal simplification allows for an en- 
hanced flexibility in the simultaneous treatment of elec- 
tric fields and strains. For example, it is possible to in- 
troduce a rigorous constant-pressure enthalpy by simply 
defining 



U*{d) = mm[U{d, n) - ntt] , 

n 



(10) 



where ir is the external pressure and fl is the cell vol- 
ume. We will demonstrate the use of this strategy in the 
application to PbTiC>3 . 

Legendre transformation. The transformation from 
variables D to variables £ can be regarded as part of 
a Legendre transformation. We spell out this connection 
here, working instead with reduced variables (^1,^2,^3) 
and (ei,£2,£3)- First, we note that 
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Recall that Si — a; • £, so that —Si is just the po- 
tential step Vi encountered while moving along lattice 
vector &i, while Qi is just the free charge on cell face 
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i. Thus, when the system undergoes a small change 
at fixed (£1,62, £3), the work done by the battery is 
— Vi dQi = — Si ddi. We therefore define 



min [U(di,d 2 ,d 3 ) 

di,d 2 ,d 3 



]eidi], (12) 



where the potentials have become the new indepen- 
dent variables and di are now implicit in the minimum 
condition. The energy functional U({di}) and ^({ei}) 
thus form a Legendre-transformation pair. 

All the gradients with respect to the internal and strain 
degrees of freedom are preserved by the Legendre trans- 
formation and need not be rederived for T . The physical 
electrical boundary conditions, however, have changed 
back to the closed- circuit case. It is therefore natural to 
expect the functional T to be closely related to the fixed- 
£ enthalpy T of Eq. ([3]). Indeed, it is straightforward to 
show that 



T = U- —£ D = T - —£ 2 

47T 8tt 



(13) 



At fixed strain and £», the Vl£ 2 /8tt term is constant, and 
thus does not contribute to the gradients with respect to 
the internal variables, consistent with Eq. ((4J. However, 
the stress derived from T differs from the one derived 
from T by the Maxwell term <r^ ax , which is absent in T 
(details of the derivation are provided in Supplementary 
Section 4.3). Although the Maxwell stress is tiny on the 
scale of typical first-principles calculations (e.g. 10 8 V/ m 
produces a pressure of 44.3 KPa), for reasons of formal 
consistency we encourage the use of T in place of J- in 
future works. 

Partial Legendre transformations. It is also possible 
to define hybrid thermodynamic functionals via partial 
Legendre transformations that act only on one or two of 
the three electrical degrees of freedom. Of most interest 
is the case of two fixed V and one fixed Q, i.e., func- 
tions of variables (e"i, £2, ^3)- The special direction is de- 
noted by unit vector q which is along direction b3 . When 
e~l — e~2 = 0, this applies to two common experimental 
situations: the case of an insulating film sandwiched in 
the q direction between parallel electrodes in open-circuit 
boundary conditions, and the case of a long-wavelength 
LO phonon of wavevector q where the q — > limit is 
taken along direction q. 

This latter case of LO phonons exemplifies the phys- 
ical interpretation of our method and its usefulness. 
While the gradients of U and its partially Legendre- 
transformed partner are identical, the force constant ma- 
trices, which are second derivatives, are not. Indeed, the 
force-constant matrices are found to differ by 



A2£ 1 a, J p = 



4tt (Zj.q) a (Z r q)j) 
O q • too • q 



(14) 



where la labels the atom / and its displacement direction 
a, Z/ Q is the corresponding dynamical charge, and eoo is 



the purely electronic dielectric tensor. This is readily 
identified as the non-analytic contribution to the LO-TO 
splitting of a phonon of small wavevector q in the theory 
of lattice dynamics Q . This demonstrates that the lattice 
dynamical-properties of a given insulating crystal within 
our fixed- D method are fully consistent with what should 
be expected from a change in the electrical boundary 
conditions from transverse to longitudinal. 

Dielectric tensor and linear response. This scheme 
lends itself naturally to the perturbative linear-response 
analysis of the second derivatives of the internal energy 
as described in Ref. with two important differences. 
First, in our scheme the derivatives at constant D become 
the elementary tensors, while the derivatives at constant 
£ are "second-level" quantities; this is an advantage, 
since using D as independent variable is very convenient 
in ferroelectric systems. Second, the use of the reduced 
field variables di and e~i in place of the macroscopic vec- 
tor fields P and £ makes the discussion of strains under 
an applied field much more rigorous and intuitive. 

As an example of the relationship between constrained- 
e and constrained-d tensors it is useful to introduce the 
inverse capacitance, 7 = C _1 , in matrix form as 



Hi 



d 2 U 
ddiddj 



(15) 



Incidentally, while this expression is fully general and 
well-defined in the non-linear regime, for the special case 
of a linear medium we can write 



(16) 



which generalizes the textbook formula U = Q 2 /2C to 
the case of three mutually coupled capacitors. It can be 
shown that the same information can be obtained within 
the constrained-e approach by means of the relationship 



d 2 T 

deidsj 



(17) 



The matrix 7^ can be thought of a "reduced" represen- 
tation of the macroscopic dielectric tensor, 



( e 1 )afJ = 7- ^2 ^ bi > a ' 
or equivalently 



£q/3 = ~n ^-^ 7 ai,a 011,13 

id 



(18) 



(19) 



We will consider, in addition to the total static ca- 
pacitance above, the closely related frozen-strain 7^ and 
frozen-ion 7?? tensors. The remainder of the response 
functions discussed in Ref. Q can be similarly defined in 
terms of the second derivatives of U ({d}, u, n). 
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FIG. 2: Potential step and internal energy as a func- 
tion of d. Symbols, calculated using constrained-D method: 
(a) reduced electric field e (squares); (b) internal energy U 
(circles). Solid curves: (a) numerical cubic spline fit to the 
symbols; (b) numerical integral of the spline fit in (a). Inset: 
enlargement near the minimum, also showing magnitude of 
the error made if Pulay stresses are neglected (dashed curve) . 



Applications. In the following we illustrate our method 
by computing the electrical equation of state of a pro- 
totypical ferroelectric material, PbTiC>3. Our calcula- 
tions are performed within the local-density approxima- 
tion [l(| (LDA) of density-functional theory using norm- 
conserving [ll[ pseudopotentials and a planewave cutoff 
of 150 Ry. The tetragonal unit cell contains one formula 
unit, and a 6 x 6 x 6 Monkhorst and Pack 



12l | grid is 



used to sample the Brillouin zone. The finite electric 
field is applied through a Wannier-based real-space tech- 
nique [1 31 ] - which con verg es quickly as a function of k- 
point mesh resolution [lj]; indeed, tests made with finer 
meshes showed no differences within numerical accuracy. 
We obtain an equilibrium lattice constant of a=3.879 A 
for cubic paraelectric PbTiC>3, in line with values previ- 
ously reported in the literature. Due to the tetragonal 
symmetry, the state of the system is fully determined 
by six parameters: the electric displacement d, the cell 
parameters a and c, and three internal coordinates de- 
scribing relative displacements along z. 

Starting from the relaxed cubic structure in zero field, 
we calculate the equilibrium state for ten evenly spaced 
values of the reduced displacement d, ranging from 
d=0.1 e to tf=1.0e (where — e is the electron charge), and 
relaxing all the structural variables at each d value. We 
set a stringent accuracy threshold of 10 -5 Ha/bohr for 
atomic forces and 10 -7 Ha/bohr 3 for stresses. First we 
check the internal consistency of the formalism by veri- 
fying that our calculated potential drop e coincides with 
the numerical derivative of U with respect to d as ex- 
pected from Eq. pip . The comparison is shown in Fig. [21 
where the discrepancies, of order 10 -6 Ha, are not even 
visible. This confirms the internal consistency of the for- 
malism and the high numerical accuracy of the calcula- 
tions. The minimum in Fig. [2] (b) [which coincides with 
the zero-crossing in (a)] at d = 0.725 e corresponds to a 
spontaneous polarization of P s = 0.78 C/m 2 . 
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FIG. 3: Dielectric properties, (a) Calculated inverse ca- 
pacitance in the free-stress (7, circles) and fixed strain , 
squares) limits. The points were obtained by extracting the 
symmetric 6x6 elementary response tensors by finite differ- 
ences (steps of ±0.001 were taken for each parameter) for each 
value of d; the continuous curve is the result of numerical dif- 
ferentiation of the splined potential, (b) Evolution of 7 for 
increasingly larger negative pressures. 



We note that this comparison is sensitive to the Pulay 
stress, even in the present case where our conservative 
choice of the plane-wave cutoff makes this error as small 
as ttp — 82 MPa. Neglecting such error corresponds to 
applying a spurious hydrostatic pressure of — 7rp, which 
leads to a discrepancy between the integrated potential 
(dashed curve in the inset of Fig. ^ and the calculated 
internal energy values. The agreement can be restored 
by plotting, instead of U (circles), the correct functional, 
Eq. (|10|) . for constant-pressure conditions (plus symbols). 
As such, this comparison constitutes a stringent test that 
all numerical issues have been properly accounted for, 
particularly in systems like PbTiC>3 that are character- 
ized by a strong piezoelectric response. 

Having verified the accuracy and consistency of our 
method, we now demonstrate its utility by analyzing 
the second derivative of the internal energy (or equiv- 
alently the first derivative of the potential) as a func- 
tion of d, which corresponds to the inverse capacitance 
7. The symbols in Fig. [3£a) show the linear-response 
values of both 7^ and 7, which are identical in the non- 
piezoelectric cubic limit. The numerical derivative of the 
splined potential of Fig. [2]^ a) accurately matches 7, again 
confirming the high numerical quality of our calculations. 
Fig. [21a) shows that the inverse capacitance is negative 
for < d < 0.395 [the zero-crossing point corresponds to 
the inflection point of the U(d) curve of Fig. Hfb), and 
to the minimum of e(d) of Fig. H(a)]. This is indicative 
of the fact that cubic PbTiC>3 is characterized by a fer- 
roelectric instability, which means that the U(d) curve 
has a negative curvature around the saddle point d = 0. 
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FIG. 4: Lattice-dynamical and piezoelectric proper- 
ties, (a) Longitudinal (filled symbols, continuous curves) and 
transverse (open symbols, dashed curves) optical modes of Ti5 
symmetry as a function of d at zero pressure. Selected TO 
frequencies for 7r=-4 GPa are shown as red star symbols, (b) 
Calculated piezoelectric coefficient as a function of external 
pressure n and bias V (thick curves). Results of the Landau- 
Devonshire model of Refs. [l6| and Qj] are shown as a thin 
dashed curve for comparison. 



We suggest, therefore, that the constrained-D inverse ca- 
pacitance at d = 0, while not accessible experimentally 
(since it corresponds to an unstable configuration of the 
crystal), is a useful indicator of the strength of the fer- 
roelectric instability. As such, it can play an important 
role in determining the critical thickness for ferroelec- 
tricity in thin pcrovskite films; in particular, a mate- 
rial with lower 7 should be ferroelectric down to smaller 
thicknesses, provided that the depolarizing effects due 
to the ferroelectric /electrode interface are equally impor- 
tant 



15| . Note that in our terminology one ferroelectric 



can be both stronger and less polar than another if it has 
a more negative 7 but a smaller \P S \- 

Piezoelectricity. Interestingly, the free-stress 7(D) 
curve in Fig. [3fa) shows a peculiar plateau for 0.8 < 
d < 1.0, indicative of a "softening" of the response of the 
crystal to the applied electric field. This behavior is sur- 
prising, as an electric field of increasing strength would 
rather be expected to drive a perovskite crystal further 
into the anharmonic regime, with a consequent progres- 
sive hardening of the overall dielectric response [17]. To 
analyze this effect, we start by observing that the evolu- 
tion of the fixed-strain 7^ dielectric response as a function 
of d is essentially featureless. This indicates that optical 
phonons alone cannot be responsible for the effect, and 
volume (and/or cell-shape) relaxations are crucially in- 
volved. To investigate the coupling between the volume 
and the dielectric response, we repeated our calculations 
within a negative hydrostatic pressure by using the mixed 
fixed-D, fixed- tt enthalpy defined in Eq. (fTT7|) . The results 
for the free-stress "f(D), plotted in Fig.[3](b), show a dra- 



matic influence of the external pressure on the dielectric 
response of the system. In particular, the plateau in 7(D) 
becomes an increasingly deeper local minimum, which 
crosses the 7 = axis for tt between —3 and —4 GPa; the 
local minimum is approximately centered in d = 0.925 
for all values of tt. 

A negative 7 indicates a structural instability, and 
structural instabilities in ferroelectric systems are usu- 
ally understood in terms of "soft" phonon modes. In 
order to see whether this picture applies here, we plotted 
in Fig. Ufa) the zone-center polar mode frequencies as a 
function of d. At zero pressure, the curves show no no- 
table irregularity, consistent with the smooth evolution 
of the fixed-strain response 7^. Remarkably, the exter- 
nal pressure has a negligible influence on such frequen- 
cies, which remain practically unchanged for the strained 
crystal at the same value of d, confirming our hypothesis 
that the effect is essentially of piezoelectric nature. 

Pursuing this idea, we combined the values of the po- 
tential drop V(d) with the equilibrium values of out-of- 
plane lattice parameter, and calculated the free-stress 
piezoelectric coefficient by numerical differentiation as 



h 



dc 
"dV 



(20) 



The results, plotted in Fig. dfb), show for zero pressure 
a clear peak at V ~ 0.2V, corresponding to a value of 
the internal field of about 450 MV/m. We identify this 
peak with the plateau in the inverse capacitance curve 
of Fig. [3]^a) . For incresingly large negative pressures, the 
piezoelectric peak becomes sharper and shifts to smaller 
values of the potential; for tt < —3 GPa (not shown) the 
piezoelectric coefficient diverges, corresponding to the 
crossover to the unstable region in Fig. [3] (b). 

These results shed light on the recent experimental 
measurements of Ref. Il8l where a remarkable anomaly 
in the piezoelectric response of PZT films in high fields 
{£ = 200 — 300 MV/m) was detected. Such an anomaly 
was rationalized in terms of a transition to a superte- 
tragonal state, which previous first-principles calcula- 
tions [lj| predicted to be stable in PbTiOa under an ap- 
plied negative hydrostatic pressure. However, the ques- 
tion remained whether the electric field might actually 
produce such a transition within the experimental range 
of applied fields. Our results demonstrate that the piezo- 
electric coupling is indeed capable of driving such a tran- 
sition, at least under isotropic free-stress boundary con- 
ditions. (In future work, we plan to extend these investi- 
gations by considering epitaxial strain clamping effects.) 
We note that our calculated value of the piezoelectric co- 
efficient of PbTi03 at zero field and pressure (h = 82.5 
pm/V) is in excellen t ag reement with previous Landau- 
Devonshire theories [id. Tl7j|. but the evolution of h for 
nonzero values of the applied potential substantially dif- 
fers [see Fig. H](b)]. For small values of the electric field, 
in particular, our ab-initio results indicate that h remains 
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roughly constant. Then h increases significantly at higher 
fields, up to a value £ ~ 450MV/m, where it starts de- 
creasing again. A monotonic decrease was predicted in- 
stead by the model of Ref. [T3- This result, therefore, has 
important implications for the tunability of the piezore- 
sponse of lead titanate crystals. 

Summary and outlook. In conclusion, we have pre- 
sented a formalism that provides full control over the 
electrical degrees of freedom in a periodic first-principles 
electronic-structure calculation. We have in mind several 
immediate applications for our method. First and fore- 
most, fixing D in ferroelectric capacitors by using the 
methods of Ref. [l3j will allow for a detailed analysis of 
the microscopic mechanisms determining the depolariz- 
ing field, both in the linear and anharmonic regime, an 
issue which is central to the development of efficient ferro- 
electric devices. Second, imposing constant- D electrical 
boundary conditions has the virtue of making the force- 
constant matrix of layered heterostructures short-ranged 
in real space. This allows one to accurately model the po- 
larization and response of complex superlattices, capaci- 
tors and interfaces in terms of the electrical properties of 
the elementary building blocks; a demonstration of this 
strategy was recently reported in Ref. [2(J Third, com- 
plex couplings between different order parameters can 
now be treated with unprecedented flexibility, opening 
new avenues in the theoretical study of magnetoelectric 
multiferroics and improper ferroelectrics. 
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